function [output] = ana_dynamics
%ANA_DYNAMICS	analyse dynamical changes in pulse pressure, heart rate and cardiac output
%	[output] = ana_dynamics
%
%
%	Inputs:         N/A
%
%
%
%	Outputs:        output
%
%
%	See also explore_data

%	References: N/A
%
%

%	Copyright 2013 MAF Pimentel
%	This program is free software: you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation, either version 3 of the License, or
%	(at your option) any later version.
%
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%
%	You should have received a copy of the GNU General Public License
%	along with this program.  If not, see <http://www.gnu.org/licenses/>.


%	$Author: MAF Pimentel$
%	$Revision: 1.0.0.0$
%	$Date: 31-May-2013 08:54:57$
%	Contact: marco.and.pimentel@gmail.com
%	Originally written on: PCWIN64

close all; clc;

DOSAVE = 0;
DOPLOT = 1;
dirout = '../outputs/dynamics/';


% Load ptstruct
load '../data/ptstruct.mat';
N = length(ptstruct);


%% Assign patients to each class of interest

clvef = get_lvefclass(ptstruct, DOPLOT);
cmort = get_mortclass(ptstruct, DOPLOT);

%% Determine Cardiac Output (index)
% ???


%% Get average time series
%
ismean = 1;
win = 4;
noverlap = 0;
for i = 1 : 4
    [clvef{i+1}.ppm, clvef{i+1}.hrm, clvef{i+1}.lam, clvef{i+1}.uom, clvef{i+1}.com] = get_meantimeseries(ptstruct,clvef{i+1}.id,win,noverlap,ismean);
end
for i = 1 : 2
    [cmort{i}.ppm, cmort{i}.hrm, cmort{i}.lam, cmort{i}.uom, ~] = get_meantimeseries(ptstruct,cmort{i}.id,win,noverlap,ismean);
end

% plot_dynamics1(clvef);
plot_dynamics2(cmort);


%% Perform analysis with the lactate and urine output

% number of hours before lactate measurement to average
win = 4*60; % in mins
for c = 1 : 2
    for p = 1 : length(cmort{c}.id)
        i = cmort{c}.id(p);
        N = size(ptstruct(i).data{1}.lactate,1);
        cmort{c}.la{p} = zeros(N,5);
        
        for n = 1 : N
            t = ptstruct(i).data{1}.lactate(n,1);
            indpp = find(t-win < ptstruct(i).data{1}.pulsep(:,1) & ...
                ptstruct(i).data{1}.pulsep(:,1) <= t);
            indhr = find(t-win < ptstruct(i).data{1}.hr(:,1) & ...
                ptstruct(i).data{1}.hr(:,1) <= t);
            induo = find(t-win < ptstruct(i).data{1}.uo(:,1) & ...
                ptstruct(i).data{1}.uo(:,1) <= t);
            
            cmort{c}.la{p}(n,1) = t;
            cmort{c}.la{p}(n,2) = median(ptstruct(i).data{1}.pulsep(indpp,2));
            cmort{c}.la{p}(n,3) = median(ptstruct(i).data{1}.hr(indhr,2));
            cmort{c}.la{p}(n,4) = ptstruct(i).data{1}.lactate(n,2);
            cmort{c}.la{p}(n,5) = cmort{c}.la{p}(n,2)*cmort{c}.la{p}(n,3);
            cmort{c}.la{p}(n,6) = median(ptstruct(i).data{1}.uo(induo,2));
        end
        
%         N = size(ptstruct(i).data{1}.uo,1);
%         cmort{c}.uo{p} = zeros(N,5);
%         for n = 1 : N
%             t = ptstruct(i).data{1}.uo(n,1);
%             indpp = find(t-win < ptstruct(i).data{1}.pulsep(:,1) & ...
%                 ptstruct(i).data{1}.pulsep(:,1) <= t);
%             indhr = find(t-win < ptstruct(i).data{1}.hr(:,1) & ...
%                 ptstruct(i).data{1}.hr(:,1) <= t);
%             cmort{c}.uo{p}(n,1) = t;
%             cmort{c}.uo{p}(n,2) = median(ptstruct(i).data{1}.pulsep(indpp,2));
%             cmort{c}.uo{p}(n,3) = median(ptstruct(i).data{1}.hr(indhr,2));
%             cmort{c}.uo{p}(n,4) = ptstruct(i).data{1}.uo(n,2);
%             cmort{c}.uo{p}(n,5) = cmort{c}.uo{p}(n,2)*cmort{c}.uo{p}(n,3);
%         end
        fprintf('%.1f Completed\n',p/length(cmort{c}.id)*100);
    end
end

% Columns order: [time pp hr la co uo]



%% Explore KDE dynamics
% select parameters to investigate -> t, x and y (time must come first)

% param = [1 4 6]; % For example [1 2 4] corresponds to [time pp la] or [time pp uo]
% 
% if length(param) ~= 3
%     fprintf('\nSelect 4 parameters that include class and time in the first two columns');
%     return;
% end

% win = 12 * 60;
% kde_dynamics(cmort{1},cmort{2},param,win,1);

end

function kde_dynamics(c1,c2,param,win,DOPLOTVIDEO)
% build KDE model for each one of the classes and build video with
% differences bettween the models from the two classes
if nargin < 5
    DOPLOTVIDEO = 1;
end

addpath('netlab/');

% create empty matrices for each class
C1 = []; C2 = [];
for i = 1 : length(c1.la)
    C1 = [C1; c1.la{i}(:,param)];
end

for i = 1 : length(c2.la)
    C2 = [C2; c2.la{i}(:,param)];
end
[r1,~] = find(isnan(C1)); [r2,~] = find(isnan(C2));
C1(unique(r1),:) = []; C2(unique(r2),:) = [];

figure; 
[f1, x1] = ksdensity([C1(:,2)]); [f2, x2] = ksdensity([C2(:,2)]); 
subplot(121); title('Feature 1'); plot(x1,f1); hold on; plot(x2,f2,'r');
legend('Class 1','Class 2')
[f1, x1] = ksdensity([C1(:,3)]); [f2, x2] = ksdensity([C2(:,3)]);
subplot(122); title('Feature 2'); plot(x2,f2,'r'); hold on; plot(x1,f1);
legend('Survivors','Non-survivors')

% Normalise features
meanval = mean([C1; C2]); meanval = meanval(2:3);
stdval = std([C1; C2]); stdval = stdval(2:3);
C1n = C1; C2n = C2;
C1n(:,2:3) = C1(:,2:3) - repmat(meanval,size(C1,1),1);
C1n(:,2:3) = C1n(:,2:3)./repmat(stdval,size(C1,1),1);
C2n(:,2:3) = C2(:,2:3) - repmat(meanval,size(C2,1),1);
C2n(:,2:3) = C2n(:,2:3)./repmat(stdval,size(C2,1),1);

% Build KDE model over different time windows and produce video
% Create 
params.na = 100; % grid number on x-axis. na and nb affects the shape of the contour plot!
params.nb = 100; % grid number on y-axis.
params.nlevels = 50;
[xgrid, xgrida, xgridb] = get_grid2D([C1n(:,2:3); C1n(:,2:3)], params.na, params.nb);

nwin = 0:win:min([prctile(c1.los,75) prctile(c2.los,75)])*24*60;

if DOPLOTVIDEO
    aviobj = avifile('example4.avi','compression','None');
    aviobj.quality = 100;
    aviobj.KeyFramePerSec = 1;
end
fig = figure;
lims = [-3 3];
for i = 2 : length(nwin)
    ind = find(nwin(i-1) <= C1n(:,1) & C1n(:,1) < nwin(i));
    data1 = C1n(ind,2:3);
    ind = find(nwin(i-1) <= C2n(:,1) & C2n(:,1) < nwin(i));
    data2 = C2n(ind,2:3);
    mix1 = Train_kde(data1);
    mix2 = Train_kde(data2);
    
    ygrid_vec1 = Out_kde(xgrid, mix1);
    ygrid1 = reshape(ygrid_vec1, params.na, params.nb); % ygrid is na by nb; 
    ygrid_vec2 = Out_kde(xgrid, mix2);
    ygrid2 = reshape(ygrid_vec2, params.na, params.nb); % ygrid is na by nb;
    
    
    subplot(221)
    plot(data2(:,1),data2(:,2),'or','MarkerFaceColor','r','MarkerEdgeColor','k');
    hold on
    plot(data1(:,1),data1(:,2),'og','MarkerFaceColor','g','MarkerEdgeColor','k');
    hold off; title(sprintf('Survivors, Day = %.2f',nwin(i)/(24*60)));
%     xlim([min(xgrida(:)) max(xgrida(:))]); ylim([min(xgridb(:)) max(xgridb(:))]);
    xlim(lims); ylim(lims); drawnow;
    
    subplot(222)
    plot(data1(:,1),data1(:,2),'og','MarkerFaceColor','g','MarkerEdgeColor','k');
    hold on
    plot(data2(:,1),data2(:,2),'or','MarkerFaceColor','r','MarkerEdgeColor','k');
    hold off; title(sprintf('Non-Survivors, Day = %.2f',nwin(i)/(24*60)));
%     xlim([min(xgrida(:)) max(xgrida(:))]); ylim([min(xgridb(:)) max(xgridb(:))]);
    xlim(lims); ylim(lims); drawnow;
    
    subplot(223)
    contourf(xgrida, xgridb, ygrid1, params.nlevels); shading flat; % Filled 2D contour plot.
%     xlim([min(xgrida(:)) max(xgrida(:))]); ylim([min(xgridb(:)) max(xgridb(:))]);
    xlim(lims); ylim(lims); drawnow;
    
    subplot(224)
    contourf(xgrida, xgridb, ygrid2, params.nlevels); shading flat; % Filled 2D contour plot.
    xlim(lims); ylim(lims); drawnow;
    
    if DOPLOTVIDEO
        F = getframe(fig);
        aviobj = addframe(aviobj,F);
    end
    
end
if DOPLOTVIDEO
    close(fig);
    aviobj = close(aviobj);
end

end

function plot_dynamics2(cmort)
% plot averaged timeseries for different variables included in cmort
colorspec = {[0 1 0],[1 0 0]};
LABELS = {'Survivors','Non-Survivors'};
PLOTMEAN = 1;

figure;
if PLOTMEAN
    for i = 1 : 2
        errorbar(cmort{i}.ppm(:,1),cmort{i}.ppm(:,2),cmort{i}.ppm(:,3),'Color',colorspec{i});
        hold on;
    end
    legend(LABELS);
else
    for i = 1 : 2
        plot(cmort{i}.ppm(:,1),cmort{i}.ppm(:,4),'-','Color',colorspec{i},'LineWidth',2);
        hold on;
        plot(cmort{i}.ppm(:,1),cmort{i}.ppm(:,5),'--','Color',colorspec{i});
        plot(cmort{i}.ppm(:,1),cmort{i}.ppm(:,6),'--','Color',colorspec{i});
    end
    legend(LABELS{1},'','',LABELS{2},'','');
end
xlabel('Time (days)'); ylabel('Pulse Pressure (mmHg)');

figure;
if PLOTMEAN
    for i = 1 : 2
        errorbar(cmort{i}.hrm(:,1),cmort{i}.hrm(:,2),cmort{i}.hrm(:,3),'Color',colorspec{i});
        hold on;
    end
    legend(LABELS);
else
    for i = 1 : 2
        plot(cmort{i}.hrm(:,1),cmort{i}.hrm(:,4),'-','Color',colorspec{i},'LineWidth',2);
        hold on;
        plot(cmort{i}.hrm(:,1),cmort{i}.hrm(:,5),'--','Color',colorspec{i});
        plot(cmort{i}.hrm(:,1),cmort{i}.hrm(:,6),'--','Color',colorspec{i});
    end
    legend(LABELS{1},'','',LABELS{2},'','');
end
xlabel('Time (days)'); ylabel('Heart Rate (bpm)');

figure;
if PLOTMEAN
    for i = 1 : 2
        errorbar(cmort{i}.lam(:,1),cmort{i}.lam(:,2),cmort{i}.lam(:,3),'Color',colorspec{i});
        hold on;
    end
    legend(LABELS);
else
    for i = 1 : 2
        plot(cmort{i}.lam(:,1),cmort{i}.lam(:,4),'-','Color',colorspec{i},'LineWidth',2);
        hold on;
        plot(cmort{i}.lam(:,1),cmort{i}.lam(:,5),'--','Color',colorspec{i});
        plot(cmort{i}.lam(:,1),cmort{i}.lam(:,6),'--','Color',colorspec{i});
    end
    legend(LABELS{1},'','',LABELS{2},'','');
end
xlabel('Time (days)'); ylabel('Lactic Acid');

figure;
if PLOTMEAN
    for i = 1 : 2
        errorbar(cmort{i}.uom(:,1),cmort{i}.uom(:,2),cmort{i}.uom(:,3),'Color',colorspec{i});
        hold on;
    end
    legend(LABELS);
else
    for i = 1 : 2
        plot(cmort{i}.uom(:,1),cmort{i}.uom(:,4),'-','Color',colorspec{i},'LineWidth',2);
        hold on;
        plot(cmort{i}.uom(:,1),cmort{i}.uom(:,5),'--','Color',colorspec{i});
        plot(cmort{i}.uom(:,1),cmort{i}.uom(:,6),'--','Color',colorspec{i});
    end
    legend(LABELS{1},'','',LABELS{2},'','');
end
xlabel('Time (days)'); ylabel('Urine Output');


end

function plot_dynamics1(clvef)
% plot averaged timeseries for different variables included in clvef
colorspec = {[1 0 0],[0 1 0],[0 0 1],[1 .5 0]};
LABELS = {'LVEF 1','LVEF 2','LVEF 3','LVEF 4'};
PLOTMEAN = 1;

figure;
if PLOTMEAN
    for i = 1 : 4
        errorbar(clvef{i+1}.ppm(:,1),clvef{i+1}.ppm(:,2),clvef{i+1}.ppm(:,3),'Color',colorspec{i});
        hold on;
    end
    legend(LABELS);
else
    for i = 1 : 4
        plot(clvef{i+1}.ppm(:,1),clvef{i+1}.ppm(:,4),'-','Color',colorspec{i},'LineWidth',2);
        hold on;
        plot(clvef{i+1}.ppm(:,1),clvef{i+1}.ppm(:,5),'--','Color',colorspec{i});
        plot(clvef{i+1}.ppm(:,1),clvef{i+1}.ppm(:,6),'--','Color',colorspec{i});
    end
    legend(LABELS{1},'','',LABELS{2},'','',LABELS{3},'','',LABELS{4},'','');
end
xlabel('Time (days)'); ylabel('Pulse Pressure (mmHg)'); xlim([0 25]);

figure;
if PLOTMEAN
    for i = 1 : 4
        errorbar(clvef{i+1}.hrm(:,1),clvef{i+1}.hrm(:,2),clvef{i+1}.hrm(:,3),'Color',colorspec{i});
        hold on;
    end
    legend(LABELS);
else
    for i = 1 : 4
        plot(clvef{i+1}.hrm(:,1),clvef{i+1}.hrm(:,4),'-','Color',colorspec{i},'LineWidth',2);
        hold on;
        plot(clvef{i+1}.hrm(:,1),clvef{i+1}.hrm(:,5),'--','Color',colorspec{i});
        plot(clvef{i+1}.hrm(:,1),clvef{i+1}.hrm(:,6),'--','Color',colorspec{i});
    end
    legend(LABELS{1},'','',LABELS{2},'','',LABELS{3},'','',LABELS{4},'','');
end
xlabel('Time (days)'); ylabel('Heart Rate (bpm)'); xlim([0 25]);

figure;
if PLOTMEAN
    for i = 1 : 4
        errorbar(clvef{i+1}.com(:,1),clvef{i+1}.com(:,2),clvef{i+1}.com(:,3),'Color',colorspec{i});
        hold on;
    end
    legend(LABELS);
else
    for i = 1 : 4
        plot(clvef{i+1}.com(:,1),clvef{i+1}.com(:,4),'-','Color',colorspec{i},'LineWidth',2);
        hold on;
        plot(clvef{i+1}.com(:,1),clvef{i+1}.com(:,5),'--','Color',colorspec{i});
        plot(clvef{i+1}.com(:,1),clvef{i+1}.com(:,6),'--','Color',colorspec{i});
    end
    legend(LABELS{1},'','',LABELS{2},'','',LABELS{3},'','',LABELS{4},'','');
end
xlabel('Time (days)'); ylabel('Cardiac Output = HR x PP'); xlim([0 25]);
end
